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' We consider a gas of N(= 6, 10, 15) Bose particles with hard-core repulsion, contained in a 

' quasi-2D harmonic trap and subjected to an overall angular velocity Q, about the z-axis. Exact 
diagonalization of the n x n many-body Hamiltonian matrix in given subspaces of the total (quan- 

^ ' tized) angular momentum L^, with n ~ 10^ (e.g. for L2=N=15, n =240782) was carried out using 

Q , Davidson's algorithm. The many-body variational ground state wavefunction, as also the corre- 

■ sponding energy and the reduced one-particle density-matrix p(r,r') = '^f X/j (r) Xm (i"') were 
calculated. With the usual identification of Q as the Lagrange multiplier associated with for a 

, rotating system, the Lz — f2 phase diagram (or the stability line) was determined that gave a num- 

' ' ber of critical angular velocities Qci, i = 1, 2, 3, • • ■ , at which the ground state angular momentum 

I and the associated condensate fraction, given by the largest eigenvalue of the reduced one-particle 

'~~] . density-matrix, undergo abrupt jumps. A number of (total) angular momentum states were found 

' to be stable at successively higher critical angular velocities fid, i = 1,2,3, •■ • for a given N. All 

, the states in the regime N > Lz > are metastable. For >N, the values for the stable ground 

■ states generally increased with the increasing critical angular velocities Qci, and the condensate was 
' strongly depleted. The critical fid values, however, decreased with increasing interaction strength 

C , as well as the particle number, and were systematically greater than the non-variational Yrast-state 

. ' values for the single vortex state with Lz=N. We have also observed that the condensate fraction 

I for the single vortex state (as also for the higher vortex states) did not change significantly even as 

a the 2-body interaction strength was varied over several(~ 4) orders of magnitude in the moderately 

I ' to the weakly interacting regime. 

^ ' PACS numbers: 05.30.Jp, 67.40.-w, 03.75.Fi 

^ ■ 

o ■ 

, CJ, . I. INTRODUCTION 

• The-|ep5perimental realization of Bose-Einstein Condensation(BEC) in dilute vapors of ultra-cold(nanokelvin) alkali 
^ ' atomgJTj, and other SYStemsQ, trapped in harmonic potential wells has qualitatively extended the domain of occurrence 

^ ^ . of the quantum fluidso 13. Unlike their dense and strongly interacting homogeneous (bulk) counterparts, e.g. liquid 
' ^He, these mesoscopic gaseous systems are dilute, weakly or moderately interacting and inhomogeneous— with conto 
, lable density, effective dimensionality, and tunable atom-atom interactions of either signQ. Further, the creationli 

^—i of the vortex sitates with quantized circulation in externally rotated traps, as also the direct observation of phase co- 
\ herence effectstll, has clearly revealed the phase rigidity characteristic of superfluidity associated with the BEC. The 

• BEC in a weakly and repulsively interacting dilute Bose gas (i.e. with the s-wave scattering length a^c <C the mean 
inter- atomic spacing and with the number of atoms in the condensate 1) has often been described macrosconi- 
cally through the Gross-Pitaevskii equation based on the condensate amplitude as a slowly varying order-parametem. 
Microscopic treatments going beyond the meanfield. axinroximation and based on the many-body variational wave- 
functions also exist, but involve heavy computationsllallj even for modest size of the system (N ~ 10 — 50). Besides, 

^ ' these many-body calculations have involved only the "Lowest Landau Level(LLL)" single-particle orbitals, with only 
Q ] the positive sign(with respect to the trap angular velocity ^Cz) of the single-particle angular momentum quantum 
O ' number m. The inter-atomic repulsion is known to qualitatively change the ground state properties of a, system of 
^ " weakly interacting neutral Bose gas (confined in a rotating harmonic trap) at T=0 in two distinct waysQ. First, it 
leads to a depletion of the condensate fraction, which is equal to 1 for an ideal(non-interacting) Bose gas. Second, 
it gives a phase rigidity, or stiffness, to the ground state many-body wavefunction. This is responsible for the super- 
fluid flow that manifests in the successive appearance of vortices with higher quantized circulations beyond certain 
critical angular velocities ilci of the rotating trap. Both these aspects are well described by the one-particle reduced 
density-matrix pi (r, r') obtained from the N-body ground-state wavefunction (ri, • • • , r-^) by partial integration, 

or tracing out of the N-1 co-ordinates from the N-body pure density-matrix ' ' ' j i"n) (^^[i • • • i r^^^ . (It is to 

be recalled in passing here that the condensate and the superfluid fractions are not the same thing. More specifically, 
e.g., for the non-rotating ground state, the condensate fraction is generally less than unity while the superfluid fraction, 
that includes the condensate as well as the above-the-condensate fraction, is exactly equal to unity. However, both 
the condensate as well as the superfluid fractions are characterized by the same quantum-mechanical phase whose 
gradient gives the superfluid velocity. Also, both vanish together above the critical temperature Tc- In fact, it is 
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the condensate that ampHfies the effect of even the weakest repulsive interaction leading to absence of single-particle 
excitation and giving rise to the phase rigidity) . 

In this work, we have studied the effect of the 2-body repulsive interaction on the condensate fraction, and on the 
the critical angular velocities(r2ci, * = 0, 1, 2, • • •) for the appearance of different vortex states for a quasi-2D Bose gas 
confined in a highly anisotropic harmonic trap (with uj^ S> which is subjected to an overall rotation U, about the 
z-axis. To this end we use a variational approach and calculate the ground state energy, the ground state wavefunc- 
tion and the associated one-particle reduced density-matrix for different values of the interaction, in given subspaces 
of the total quantized angular momentum L^. The many-body variational ground state wavefunction is obtained 
through the exact diagonalization of the nxn many-body Hamiltonian matrix with, e.g., n= 240782 for L^^N^IS, 
using Davidson's algorithm. A distinctive feature of the present work is that in constructing the many-body vari- 
ational ground state wavefunction for a given (= the eigenvalue corresponding to the z-component of the total 
angular momentum operator L^), we have included, in the configuration interaction, the onc-particle states with the 
single-particle angular momentum ^= i^^ quantum number m of either sign, as also the higher "Landau-Level(LL)" 
states. One of our main results is the variation of the critical angular velocities Oci with the interaction strength 
A2 = (N-1) y^^(the meanings of various symbols will be given in the next section), namely, that not only does 
VLci decrease with increasing A2, it also stays systematically higher thanjits non-varional value, e.g., that given by the 
relation fici = ijJA. (l — ^) for the weakly interacting, dilute caseoEjIlj. We have also observed that the condensate 
fraction for the single vortex state with Lz=N(as also for the higher vortex states) does not change significantly even if 
the 2-body scattering length a^c is changed over ^ 4 orders of magnitude in the moderately to the weakly interacting 
regime, namely, from a^c = lOOOao to lap, where ag is the Bohr radius. 

This paper is organized as follows. In Section 2, we begin with a more general system to bring out, in pass- 
ing, the mathematical analogy between the system under study here and the other well known systems like in the 
Landau-Darwin-Fock problem and argue that in a certain limiting case of interest to us here, it becomes essential 
to go beyond the Lowest Landau Level(LLL) approximation so as to include higher LLs and with the single-particle 
angular momentum eigenvalue m taking both positive as well as negative values in constructing the many-body ba- 
sis functions. Section 3 describes briefly the constrction of the many-body basis functions and the determination 
of the variational ground state wavefunction by diagonalizing the many-body Hamiltonian matrix using Davidson's 
algorithm. In Section 4, we outline the procedure for determining the critical angular velocities for the entry of the 
vortices into the system, and the determination of various density profiles obtained from the one-particle reduced 
density-matrix characterizing the vortical state. Finally, in Section 5, we present and discuss our results, and end the 
paper with a brief conclusion. 



II. THE SYSTEM AND THE HAMILTONIAN 



We begin by considering the general case of a system of interacting, spinless, charged particles (bosons) confined 
in an external harmonic potential(trap). The 2-body interaction potential is, however, assumed to be gaussian in 
the particle-particle separation. The trap is also subjected to an externally impressed rotation at an angular velocity 
rj = riz, and to a uniform magnetic field B = Bz. The Hamiltonian for the system in a frame co-rotating with the 
angular velocity f2 is then 
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is the electromagnetic vector potential in the symmetric 



Here q is the charge of the particles and A(r) = ^ {xcy — ye 

gauge. The co-ordinates {vi\ and the corresponding canonical momenta {r^i} refer to the laboratory frame. From 
now onwards, we will drop the superscript lab on the Hamiltonian H , the total angular momentum L^ , the single- 
particle angular momentum 1^ and their respective eigenvalues, and these will always be assumed to refer to the 
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laboratory frame unless otherwise specified. 

We now make the following assumptions. The confining asymmetric harmonic potential is highly oblate spheroidal, 
with Xz = ^ ^ 1, and hence our confined system is effectively quasi-2D with the x-y rotational symmetry; the 



repulsive 2-body scattering is dominantly in the s-wave channel with a scattering length Usc and 



having the 



dimension of Energy x Volume. The range a of the 2-body interaction is small enough compared to the inter-atomic 
spacing so as to effectively give a 5-function interaction potential V (r, r') ~ 4^Kn_a^g _ ^ly 

Now, for our N-body variational calculation we need to construct the N-body basis-functions with proper symmetry. 
Since the system above is rotationally invariant in the x-y plane, the z-component of the total angular momentum(L2) 
is a good quantum number leading to block diagonalization of the Hamiltonian matrix into the subspaces of L^. 
The N-body basis-functions are, in turn, constructed as linear combinations of the symmetrized products of a finite 
number of single-particle basis-functions, which are chosen to be eigenfunctions of the unperturbed single-particle 
Hamiltonian. With this in mind, let us consider the non-interacting single-particle Hamiltonian h in the rotating 
frame, which we now separate into the z and the (x-y)_plane (commuting) components, and h_L respectively: 
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The eigensolutions for hz are then given by 

^z^n. {z) = e„,u„^ (z) , where e„^ 



v^2"-nz 



e-5«»- H„, {a^z), 



and = 



nz + 2 I T^z, riz = 0,1,2,- 



= inverse of the longitudinal oscillator length (a^) 



Here H„^ is the Hermite polynomial. We assume the system to be quasi-2D in that there is practically no excitation 

along the relatively stiffer z-axis, and hence we set = 0. 

Let us now consider h_L. Using ft — QBz, B = Be^ and rj^ — xe^ + yey, we re- write it more explicitly as 

-^"^'18/ d \ 1 ^2 "1 
_r_i_dr± \'^dr±) r']_d(jp'_ 
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I +Lo'^, , the frequency for the harmonic confinement in the x-y .plane, 

2m / 

- 1 , the cyclotron-rotational(or the centrifugal) angular velocity 



about the z-axis. 

Here IC± is the single-particle 2D-harmonic oscillator Hamiltonian, for which the eigensolutions are known to be: 



with u„,„(rj_,0) 
where a± 



ik[n-\m\]y. 
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the inverse of the transverse oscillator length (aj,) , 
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2nr+ \ m\ +l\nC,^ = {n+l) hC±, 
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with rir = 0, 1, 2, • • • , m = — cx), •••,—1,0, +1, • • • , +cx), 
or equivalently, n — Q, 1, 2, •••, ni — +n, +n — 2, • • • , — n + 2, — n. (4) 



Here Li^|j_|^|-, (Q!'!''!) is the Associated Laguerre polynomial. It may be noted in passing that the above Hamiltonian 



represents the well known Landau-Darwin- Fockll3ilZl problem in which is the frequency for the harmonic confining 
potential in the x-y_plane and is the cyclotron- rotational (centrifugal) angular velocity about the z-axis. Limiting 
to — and taking m = 0, +1, +2, +3, • • • , corresponds to the "Lowest-Landau Level(LLL)" approximation. 
Let us examine the centrifugal/mechanical stability of the above system by re- writing hx as 

h± = ^ f Y^_L - (B X r±) - m (fi x r^) 



\^m.<J'^T\ - (f7 X r^)^ - (B x r_L) • (fJ x r^) = f + U 



Here T is clearly a positive-definite operator. It can readily be shown that hi is negative-definite, null or positive 
definite according as (C^ + <;) ~ is negative , zero or positive, respectively. 

For lA negative-definite, the Hamiltonian hj^ = T + is unbounded from below and there are no stable solutions. 
In the simpler case of vanishing magnetic filed B=0, this situation arises for wj^ < f2, i.e, when the rotational angular 
velocity f2 becomes larger than the confining harmonic trap frequency wj^. 

For the special case of U null, we have C± = ^ and the Hamiltonian reduces to hj^ = /C^ — C± U • This gives 
rise to a situation analogous to the Landau problem when the frequency for the harmonic confinement is equal to the 
centrifugal frequency. Setting C,±^ = 5^c, we have the eigenvalue equation 

with (e„,m - mhC±) 



( \ 

+ - (I m I -to) +- 



hcjc = ( TV + i ) hcjc 



V :v J 

where rir = 0, 1, 2, 3, • • • , and to = — oo • • • — 2, —1, 0, +1, +2, • • • + oo; 
or equivalently A/" = 0, 1, 2, 3, • • • , 

and TO = -A/", -(TV -!),•• •,-2,-1, 0,-^1,-^2, •••,-!-{». (5) 

Each of the TV levels, the so-called Landau Levels(LL), is infinitely degenerate corresponding to the infinitely many 
possible values of to. It is clear from the ordering of the single-particle energy levels (without interaction) in equa- 
tion (^ that the single-particle states with positive to values (i.e those with the angular momentum parallel to the 
overall rotational angular velocity ^ = (^i-z) are energetically favored. These states constitute a massively degenerate 
manifold. This is the usual rationale for using the positive to values only, in constructing the variational wavefunction 
for the Landau-like problem. This degeneracy for the special case oi (± ~ is, however, lifted by the inter-particle 
interactions. 

Finally, we consider the physically interesting case oiU positive-definite. The single-particle non- interacting Hamil- 
tonian hj_, now becomes hj^ ICj_ — ^ Iz with the eigenvalue solution: 

hj^u„^,„ (f±» = (^_L - dz) u„^„ (r_L, (f)) = (e„,™ - mh<;) u„^„ (rj_, 0) 

(e«,m - rnM) = 2 + ^ ^| to | -m-^^ + fiC±, 

where = 0, 1, 2, • • • , and to — — cxd, •••-2,-1,0, +1, +2, • • • + oo. (6) 

In this physically relevant situation, the particle as observed in the rotating frame finds itself in a shallower harmonic 
potential of frequency \/ C± ~ a-nd the states will be more spread out. It can be seen from equation(H) that for the 
centrifugal frequency <, significantly smaller than the confining frequency C_l , the degeneracy of the Landau levels is 
lifted even without interaction. Further, the interaction between the particles causes the the different single-particle 
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angular momentum states to scatter into each other. Thus, for the slowly rotating systems and for moderately, and 
a fortiori for strongly interacting systems, it may be necessary to include single-particle basis functions with different 
values of rir and with the angular momentum quantum number m taking both positive and negative values. This is 
the case we shall be concerned with in the following with B=0 and q = 0. 

III. THE MANY-BODY VARIATIONAL WAVEFUNCTION 

We need to perform the diagonalization of the Hamiltonian matrix in the subspaces of only. The N-body 
variational wavefunction \J/ (ri, • • • , rjyf) is a linear combination of the symmetrized products tpb (ri, ^2, ': rjvg-) of 
the one-particle basis functions Un, m, Uz (r) = Un (r) introduced earlier: 

n 

*(ri, r2,---,rN) =^ Cbi'biri, r2,---,r]vf) 
b 

^ / 1 1 '^o + vi 

with r2,---,rN) =^1]P ^n"o(^*)'7^ 11 «i (^O ' • • 




j=i/o+i'iH hi^k-i + l 
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Here {Cb} are the variational parameters, and P permutes the N particle co-ordinates. Also, Vj is the occupancy of 
the jth single-particle basis function(uj). The many-body quantum index 6, labelling the many-body basis function 
V'6 (ri, • • • , rjyf) , stands for a set of single-particle quantum numbers required to satisfy the above two defining relations 
between the single-particle occupation quantum numbers {fj}, the single-particle angular momentum quantum num- 
bers {rrij}, the number of particles N and the total angular momentum Lj,. At this point it becomes more convenient 
to switch over to second-quantization notation in the occupation number (z/j) representation. The basis function is 
then 



k k 



(4) («l) '(4) "- (<) ' 
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{i^o • • • i/k) : ^ J<j = N, ^ t-jmj = \ , (7) 
j=o j=o / 

and the many-body Hamiltonian written in the second- quantized notation is 

I* = (i I 1 I j) + \Y1 ^^1''2 I 2 I ji,j2) (ala^.ala^^ - <5i,jia?^aj,) 

i, j ii, ji 12, h 

Here aj (aj) are the usual bosonic creation(annihilation) operators, and (i | 1 | j) and (ii, i2 | 2 | ji,j2) are the one-body 
and the two-body matrix-elements respectively over the single-particle basis functions. The evaluation of matrix- 
elements for the kinetic energy and the harmonic trapping potential over the single-particle basis chosen in the 
previous section is trivial. It is also possible to obtain closed form expressions for the matrix-elements for the two- 
body (as also for the n-body, n=2,3,- • •) gaussian potential(s) over the above basis as they reduce to multi-dimensional 
gaussian integrals. (In the calculations presented here, however, we have not considered the three or the higher-body 
potentials). 

The variational parameters {C"} for the ground state wavefunction ^'o = Y^b^b 'Pb ^^'-'^ determined by min- 
imizing Ko {^'o} = (^'o |H| 4*0) — Eo (*o I ^'o) with respect to ^fo- The Lagrange multiplier Eq is identified as the 
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variational energy for the ground state. The zth excited state = J^b ^I'^'b is determined by carrying out the varia- 
tional minimization in the restricted Hilbert subspace which is orthogonal to the (i — 1) states determined earlier, i.e. 
by minimizing K, {^',1 = (*,|H|*,> -E,(*, I*,) pdth | = for j = 0, 1, •■•,(*- !)• 

The Davidson algorithm of iterative diagonalizationEJ is based on a procedure where one keeps a minimum set of i 
orthogonal, trial wavefunctions which spans a small subspace iSj of the full many-body Hilbert space: 

Si EE l^-j I *j = X! ^b'^'' ' J = 0, 1, • • • * and (*j I ^k) = Sjk for j. A: = 0, 1, • • • , i | . 

This small subspace Si is chosen to contain dominant contributions from the ground state and the first few excited 
states. An i x i representation i?^*' of the Hamiltonian H is obtained over this small subspace to set up the small 
eigenvalue equation H^'^^a.i, = A^a^/, = 0, 1, 2, • • • i. iJ^'^ is an effective Hamiltonian for the system over the small 
subspace Si spanned by {'^j, j = 0, 1, 2, • • • , i}. The eigenvalue AJ, is the i^th approximate eigenvalue. The convergence 
for the vth state is achieved when the residual vector for the z/th state | A^t^) — H | '^^) — AJ, | 'i/^) becomes a null 
vector. If the convergence has not been achieved, the residual vector A^*^, after orthonormalization, is added to 
the list of trial vectors to augment the subspace Si to obtain iS^+i = {vfj, j = 0, 1, 2, • • • , i -I- 1}. The procedure is 
continued till the convergence is obtained. In the process, if the size of the subspace Si becomes unwieldingly large, 
a certain number of higher eigenvectors are dropped from Si and the procedure once again initiated. 

For N=15 particles, for example, we have carried out calculations for all the total angular momentum states in 
the regime < Lz < 3N. Diagonalization of the Hamiltonian matrix is performed for each of the subspaces of Lz 
separately. (We have set = in the single-particle basis function u„, (z) as discussed earlier.) The single-particle 
basis Un,m {ri^, (j)) with n = 2nr+ \ m \ and = 0, 1, • • • , & | m |= 0, 1, 2, • • • , spanning the 2D x-y plane for a given 

subspace of L^, is chosen as follows. It is convenient to define / = ^ , where for real x, the symbol [x] denotes the 
greatest integer less than or equal to x. For very weakly interacting particles, almost all the particles will condense 
into a single one-particle state, n = m = Z, a Yrast-like state. As the interaction becomes stronger, the particles 
start scattering out to other single-particle states around this state, i.e. some of the particles go to states with higher 
angular momentum while some of them go to states with lower angular momentum(the two-body interaction conserves 
the total angular momentum). The single-particle angular momentum for the basis functions is now chosen to be: 
m = I — rib, I — Ub + I, ■ ■ ■ I + Jib — I, I + rib, where rib is some positive integer that we have chosen to be 3, 4 or more 
depending on the strength of the interaction and the computational resources available (nf, is a kind of the size of the 
single-particle basis chosen for the calculation for a given value of and describes configurational interaction). In 
all our calculations presented here we have taken = 0, 1, and rib — 3. Thus for example, for N=15 and for the 



chosen subspace Lz = 33, we get I 



= 2, and the single-particle angular momentum quantum number takes 



values m — —1, 0, +2, +3, +4, +5. Then, with rir = 0, 1, the single-particle basis set turns out to be 

{^^0,0, Ul, + 1, "2, +2, U3,+3, U4,+4, U5,+5, "l-l, ""2,0, "S. + l, U4,+2, M5,+3, "3 -l} • 

Thus, the N(=15)-body basis functions {ijjb}, for the ~ 33 subspace, are to be constructed from these single- 
particle basis functions. These were used in turn to construct the variational trial function = J^b ^b4'b, and hence 
the ground state properties for the given value of . 



IV. CRITICAL ANGULAR VELOCITIES AND DENSITY PROFILES FOR THE CONFINED ROTATING 

BOSE GAS AT ZERO TEMPERATURE 

From the variational ground state wavefunction (ri, • • • ,rp^) obtained in the previous section, we now, in the 
following, go on to calculate various quantities of interest. 

Let us first calculate the critical angular velocities {r^ci} for the onset of different vortical states. Note that for the 
system of N particles confined in a trap rotating at an angular velocity Q, the thermodynamic equilibrium corresponds 

to the minimum of the free energy F given by exp {— /3F} = Tr ^exp (^H — hflh^ ) }) ' ^-'^^'^^ ^ ~ ^^^z = 
as we have noted earlier, is the Hamiltonian of the system in the co-rotating frame. For a system at T=0, the 

\ /^--Ja-bX / - lab 



expression for the free energy reduces to a more simpler form F = ^^Pq ^ ^ j \ / \ 

where is the ground state wavefunction of the system in the laboratory frame for a given value of Lz obtained. 



in our case, variationally. The above relation can as well be seen as the minimization of the energy ( H ) — Eq*^ 



■ lab \ 
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subject to the constraint that the system has the angular momentum expectation value ^L^ ^ The angular velocity 

n is then the corresponding Lagrange multiplier. Since we have constructed our variational wavefunction to be 

an eigenfunction of the total angular momentum operator , we will necessarily have (Li^ ^ = Lj;. Initially, when 

the system is subjected to no external rotation, the angular momentum state = corresponds to the ground state 
of the system. As the system is rotated, other non-zero angular momentum states (L^ ^ 0) successively become the 
ground state of the system. These are obtained by minimizing the energy E™* (L^, Q, A2) in the rotating frame: 

ES°* (L„ A2) = (L„ A2) - hnh'f (8) 

Here A2 measures the two-body interaction strength and has been defined in an earlier section. Thus, the critical 
angular velocity flc beyond which the higher angular momentum state , say, becomes lower in energy in the rotating 
frame compared to the lower angular momentum state (< L^) is given by: 

_ E^CLz, A2)-Er(L^, A2) 

This gives us a number of critical angular velocities {fid, i = 0, 1, 2, • • •} at which the ground state of the rotating 
system changes its angular momentum state abruptly. The angular momentum L2 of the ground state vs the angular 
velocity fl relation, called the — f2-phase diagram, or the stability line, for the rotating system was calculated in 
terms of the variational ground states obtained. We present this in Fig. 1. 

Next, we calculate the density profile and the circulation for the vortex states. For this we need the single-particle 
reduced density-matrix. From the many-body ground state wavefunction ^'o (ri , • • • , r-^) , obtained through the 
variational exact diagonalization(ED), the one-particle reduced density matrix pi (r, r') is obtained by integrating out 
the N-1 co-ordinates: 

pi (r, r') = J J ■ ■ ■ J i"3 ■ • • i-N*o (i"' i"2, ra, • • • , r^) ^'o (r', r2, rg, • • • , r^^) 

~ ^ ^ ^ ] ^ ] ^ ] ^ ] ^ ] Pnmn^,n'm'n'^ 'Un,m,n^ i^) "^n' ,m' ,n'^ (l" ) 
n m n' m' rtj. 

= H X*c (r) (i"') , 

with X/.(r) = W ^^'^ JI'^A'^l' (10) 

n rn ^ 

where {A^j} are the eigenvalues and {x^ (r)} the corresponding eigenvectors of the one-particle reduced density-matrix 
pi (r,r'). The diagonal part pi (r,r) = p(r) gives the single-particle density-profile of the system(the condensate -I- 
the non-rotating fraction + the above-the-condensate fraction). For our finite system, the BEC corresponds to a 
single eigenvalue being significantly larger than the rest of the eigenvalues. As noted before, our system is quasi- 
2D, i.e., there is practically no excitation in the z-direction. So we take n^ = and the summation over goes 
off. Also, we note that (r) is an eigenvector of the single-particle angular momentum 1^ with eigenvalue to^, 
hence the summation over m in equation(^0|) too goes away. The only summation that we are left with is then over 
(= 2nr^ + |m^| , where n^^ = 0, 1, 2, • • •). Thus we get 

Xa'W = 4i^,m^fi^n^,m^.,o{v) (11) 

Note that n^ increases in steps of 2, which comes from quantization. Tracing out the z co-ordinate, the density 
p{rj_,(j)), the current j,p{r±,4'), the velocity v^{r±,(j)) and the circulation K{r±) profiles in the x-y plane, for the 
quasi-2D system, turn out to be 

p(r^,^)=4-^""'"^^E V 
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and j^^ (r^,0) = (r^» = 0; 



huj 



V0 [r±), 



and v^^ (r_L,(/)) = (r^, 
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V • = / V0 (rx) r^d'/' = 27r r_L V4, {r±) 

r±— constant J Q 

2 



2tt— rUfj, Ap 

TO '^-^ 



These are plotted in Figs. 2-4. For completeness we also give the velocity curl 



V X V = — ^ Ap 



p{r±,(j)) \r [ a [ fiuj 

Mn^,m^,0 - \J {nf,- I TO^ |) | TO^ |) 



X 



C.C 



(12) 
(13) 

(14) 



(15) 



(16) 



The expression for circulation as in equation(|15|) calls for some discussion. The circulation associated with each of 
the eigencomponents(x^) of the one-particle density-matrix is, of course, quantized as integral multiples of ^ for any 
closed path, and clearly can have no r± dependence for the circular paths chosen for our line integral. The expression 
in equation([l5|), however, involves an average over different components x^, and can, therefore, assume non-integral 
values that vary with r± for the chosen circular path r^=constant. This reflects the radial variation of the relative 
fractional weights A^ \Xfj. ir±, (/))|^ of the different fractions with r±. 



V. RESULTS AND DISCUSSION 

The results and discussions given below refer to the following choice of parameters: aj_ = 1.222/iTO, A^ = -s/S, 
Use = lOOOao, where ao is the Bohr radius. Please note that the scattering length a^c chosen above is 10 times more 
than as given intj for ^^Rb. Further, we present here the results for N=15 only in the angular momentum regime 

< L2 < 3N. We have, however, done the calculations for N=6 and N=10 also. First we note that we are always in 

(\ 3 1 _ i 

j <C 1, with Qho = (aj_az) ^ — a±Xz " . Further, we are also in the regime 

of moderate (a^c = lOOOao) to very weak(asc = lao) interaction strength as the parameter is of order 1 or much 

much less than 1 for the two regimes respectively. Moreover for Osc = lOOOao, the healing length ^(= {SiTnasc) ^, 

with n ^ the mean particle density) < the oscillator length Uho ^ the size of the system. Thus our study for the 

above choice of paremeters may be relevant to bulk systems. 

The critical angular velocity Oci for the single-vortex state(Lz = N) decreases for increasing value of N. Thus for 
Use = lOOOao, it is 0.88487a;^ and 0.84493w_l for N=10 and 15 respectively. For a given N, the critical angular velocity 
flci is also found to decrease with inceasing scattering length asc- Thus for N=15, r2ci(in units of uj^) are for found 
to be .99978, .99783, 0.97901 and 0.84493 for Usc = lao, lOao, lOOao, lOOOao respectively. The condensate fraction 
for a given N and corresponding to one of the stable vortex states, however, seems to be quite insesitive to a change 
in the scattering length. Thus for L^^N^IS state and age — lao, lOao, lOOao, lOOOao, the condensate fractions(Ai) 
are found to be 0.87826767, 0.87822606, 0.87777931, 0.87069438 respectively. 

In Table 1, we summarize our results for the ground state of the N=15 system in different subspaces of the total 
angular momentum in the regime < < 3N for age = lOOOao. As the system is rotated, a number of (total) 
angular momentum states are found to be stable at successively higher angular velocities ild, i = 1, 2, • • •. We have 
also given the density-matrix eigenvalues! A^} and the corresponding single-particle angular momentum quantum 
numberjTO^} for the largest three fractions(Ai > A2 > A3). We note that, corresponding to each of the stable ground 
states of the rotating system, namely, — 15, 26, 30, 33, 35, 36, 45, • • • for N=15, the second largest fraction is 
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found to be non-rotating, i.e. m2 = 0. The single-particle angular momentum quantum number mi, corresponding 
to the largest fraction, is taken to be the vorticity of the condensate. 

Figure 1 gives the — f2 phase diagram, or the stability line. One can readily identify the successive critical angular 
velocities {rici, * = 0, li 2, • • •} for the transition from one stable state to another stable state. Several points are to be 
noted here. The non-integral values for the angular momentum per particle clearly indicates the incomplete condensa- 
tion, i.e., the fact that more than one eigenvalues {A^} of the reduced density-matrix are non-zero. This is, however, 
not to be taken as the fragmented condensateE3, but merely as depletion of the condensate/superfluid fraction due to 
the interaction and the rotation. Also, we have observed that the critical angular velocity Oci decreases monotonically 
with increasing repulsive interaction and the particle number. This can be readily understood physically from the 
expression in equation(^) for the critical angular velocity (L^) if one remembers that the non-rotating (L^ = 
angular momentum) state is more compact and, therefore, has its energy raised (because of the repulsive interaction) 
relative to the expanded higher angular momentum states. More importantly, however, our critical angular Kclpcity 
Vtci for a given N, is systematically greater than the non-variational value based on the Yrast-like stateoE3ll3'E3. We 
expect this to be due to our more accurate and, therefore, lower variational ground state energies. 

In Fig. (2a) we have plotted the total (the condensate-|-the non- rotating component -l-the above-the-condensate 
fraction) density profile p{rx_) for the L^=N single- vortex angular momentum state of the system. We have also 
plotted separately the density profiles for the condensate fraction Ai |xi (r)| corresponding to the largest eigenvalue 
Ai {— 0.87069) in Fig. (2b), the non-rotating component A2 \x2 (r)|^ corresponding to the second largest eigenvalue 
A2 (= 0.06576) in Fig. (2c), and the above-the condensate fraction in Fig. (2d) corresponding to the remaining compo- 
nents of the one-particle reduced density matrix: 



p(r,r') = Aixt(r) Xi(r') + A2x;(r)x2(r') + ^^^xl{v) X,{v') 

l-t—3 



the condensate the non-rotating com 
with mi = 1 -ponent with TO2 = the above-the 

-condensate 

where, as we have earlier, m^ is the single-particle angular momentum quantum number associated with the /ith 
component of the density-matrix. As we go to higher Q, the condensate depletes— with the non-rotating component 
and the above-the-condensate fractions becoming more pronounced(as can also be seen from Table 1). 

Figure 3 depicts the velocity profile for the total system(3a), for the condensate(3b), for the non- rotating compo- 
nent (3c), and for the above-the-condensate fractions (3d) respectively. 

Figure 4 depicts the circulation profile for circular paths of radii r± in the x-y plane: (4a) for the total system, (4b) 
for the condensate fraction, (4c) for the non-rotating fraction, and (4d) for the above-the-condensate fraction respec- 
tively. As can readily be seen, the condensate fraction has an integral circulation and is constant over space, while the 
total as well as the above-the-condensate fraction has non-integral circulation and varies with r±_ as explained earlier 



below equation(15) 



In conclusion, our many-body variational calculation for an admittedly finite number (N < 15) of particles explicitly 
resolves the structure of the rotating BEG at T=0 in terms of the various components of the reduced one-particle 
density-matrix. Further work is in progress to consider spin-1 Bose particles, where we expect qualitative changes 
in the structure of the rotating BEG. In particular, we expect a dramatic increase in the critical angular velocities 
because of the spin-polarization drag effect inasmuch as the angular momentum may be absorbed preferentially in the 
spin rather than the orbit|al| |degrees of freedom without costing kinetic energy. This could lead to spin polarization 
due to rotation of the trapE2l'E3. We expect this effect to be more pronounced here than in a classical rotating system. 
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Table caption 

Tabic 1. Gives the largest three condensate fractions Ai > A2 > A3 and the corresponding single-particle angular 
momentum eigenvalues mi , m2 and 7713 in the one-particle reduced density-matrix for the ground state of the rotating 
BEG, in given subspaces of total angular momentum L^. Also given are the ground state energies Eq''' (L^) in the 
laboratory frame and the critical angular velocities flci corresponding to the stable ground states of the rotating 
system in the angular momentum regime < < 3N. 

Figure captions 

FIG.l. Gives the plot (solid line) for the total angular momentum in units of h versus the critical rotational velocity 
f2ci(iii units of a;j_) for the rotating BEG. Dashed line is for the non-interacting case, included here for reference. 

FIG. 2. Depicts the normalized particle-density in units of aj^ for, (a) the total system, (b) the condensate fraction, 
(c) the non-rotating component, and (d) the above-the-condensate fraction, in the x-y plane. Distances are in units 
of a±. 



FIG. 3. Depicts the azimuthal velocity profile {r±) in units of y for, (a) the total system, (b) the condensate 

fraction, (c) the non-rotating component, and (d) the above-the-condensate fraction, in the x-y plane. Distances are 
in units of a±. 




FIG. 4. Depicts the circulation K{r±) = / v • dl profile in units of ^ along circular paths of radius r± in the x-y 
plane for, (a) the total system, (b) the condensate fraction, (c) the non-rotating component, and (d) the above-the- 
condensate fraction. Distances are in units of a±. 
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